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Abstract 

The discovery of charged Higgs bosons is of particular importance, since their ex- 
istence is predicted by super symmetry and they are absent in the Standard Model 
(SM). If the charged Higgs bosons are too heavy to be produced in pairs at fu- 
ture linear colliders, single production associated with a top and a bottom quark is 
enhanced in parts of the parameter space. We present the next-to-leading-order cal- 
culation in supersymmetric QCD within the minimal supersymmetric SM (MSSM), 
completing a previous calculation of the SM-QCD corrections. In addition to the 
usual approach to perform the loop integration analytically, we apply a numerical 
approach based on the Bernstein-Tkachov theorem. In this framework, we avoid 
some of the generic problems connected with the analytical method. 



1 Introduction 



The discovery of charged Higgs bosons (-ff ± ) would be instant evidence for physics beyond 
the SM, which only accommodates a single neutral Higgs boson. Charged Higgs bosons 
appear in models with two Higgs doublets, as are required by supersymmetric (SUSY) 
extensions of the SM. For this reason, there is much interest in charged-Higgs-boson 
physics (review articles on this subject include, for instance, Refs. [IJEIE]). 

Phenomenologically, the Large Hadron Collider (LHC) will be the first collider with the 
potential to discover the H ± bosons. If the charged-Higgs-boson mass m% is not too large, 
i.e. if rri^j < mt — mb, the dominant production channel is via top-quark pair production 
gg —> tt with subsequent decay of a top quark or antiquark into a charged Higgs boson, 
t — > bH + or t — > bH~, respectively. If the H ± bosons are too heavy for this subsequent 
top-quark decay, then the dominant production channel would be bottom-gluon fusion, 
gb — > tH~ and gb — > iH + [H El El [7J [8] . In all these processes, the charged-Higgs-boson 
signal has to be carefully separated from large SM-QCD background at hadron colliders. 

Despite the fact that charged Higgs bosons may be discovered at the LHC, a precise 
determination of their properties will only be possible at a linear collider, such as the 
proposed international linear collider (ILC). If the value of m H ± is not too large, i.e. if 
vtlh± < \fs/2, the dominant production channel will be e + e~ — > H + H~ [9], [TD]. On the 
other hand, this production channel may not be accessible at the ILC because of the 
limited center-of-mass energy a/s. In this case, charged Higgs bosons may be copiously 
produced singly via the two channels [3J 



The first production channel ([ID was investigated in Ref. [3J, where single charged-Higgs- 
boson production processes were systematically compared with each other with leading- 
order (LO) accuracy. Here, we focus on the second production channel ([2]), which is 
enhanced in parts of the parameter space. Since QCD corrections are typically large, we 
present a computation with next-to-leading-order (NLO) accuracy, to 0(a s ). 

One part of the NLO calculation consists of the SM-QCD contribution, i.e. the purely 
gluonic corrections, which were already presented in Ref. Here, we add the SUSY- 
QCD contribution, due to squark and gluino loops, in the minimal SUSY SM (MSSM) [121 
[T3l [H]. Since the MSSM makes no assumption about the SUSY-breaking mechanism, 
but just uses explicit SUSY-breaking terms in the Lagrangian, it may be considered as 
representative for a wide class of SUSY models. The SM-QCD and SUSY-QCD parts 
taken together yield a complete prediction with 0(a s ) accuracy. 

In the calculation of the virtual corrections, we encounter loop integrals. Using the 
conventional analytic approach of Ref. [15], all scalar loop integrals can be expressed 
in terms of logarithms and dilogarithms. Furthermore, using the reduction algorithm of 
Ref. [16], all tensor loop integrals, i.e. integrals containing loop momenta in the numerator, 
can be expressed in terms of scalar integrals. Therefore, a full analytic solution for one- 
loop integrals exists. However, in general, this approach has a number of drawbacks. First 
of all, the number of dilogarithms in the analytic expression of a scalar integral increases 
rapidly with an increasing number of external legs. This may lead to cancellations for 
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e + e" -> tbH~ + c.c. . 
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Figure V. Feynman diagrams of the process e + e — > tbH at Born level. Those of the 
charge- conjugated process e + e~ —> btH + are not shown. 

multileg integrals in certain kinematic regions [17] . Furthermore the tensor reduction 
of Ref. [TB] introduces inverse Gram determinants. These may vanish at the phase- 
space boundary, even though the tensor coefficients themselves remain regular there. 
There may thus be cancellations among terms in the numerator, which may lead to 
numerical instabilities. While these generic problems may still be overcome in the case 
under consideration here, where at most box diagrams occur, by exercising care in the 
phase-space integrations, they become quite severe for five and more external legs. To 
address these problems, several improved reduction algorithms have been constructed [TBI 
fT9[ l20l I2T] allowing a numerically stable evaluation of the tensor integral coefficients. A 
different solution, typically used in the context of nondiagrammatic methods to calculate 
loop amplitudes, is to employ high-precision arithmetics for potentially unstable phase- 
space points (see Refs. [221 [231 [21] for reviews of these techniques). However, these 
improvements come at the price of either increased complexity of the more elaborate 
reduction algorithms or increased runtime from the high-precision evaluations. 

Finally, within dimensional regularization in D = 4 — 2e space-time dimensions, the 
evaluation of the loop-by-loop contribution to a two-loop correction makes it necessary to 
expand the one-loop integrals beyond the constant term in the expansion about D = 4 
dimensions. An analytic calculation of these higher-order terms is rather complicated. 

We explore here an alternative numerical approach to the evaluation of one-loop tensor 
and scalar integrals. This strategy is described in detail in Ref. [25] and is based on the 
Bernstein-Tkachov (BT) theorem [26], which can be used to rewrite one- loop integrals in 
Feynman-parametric representation in a form better suited for numerical evaluation. This 
technique yields a fast and reliable numerical calculation of multileg one-loop integrals. 
This method is numerically stable also for exceptional momentum configurations and 
easily allows for the introduction of complex masses and the calculation of higher orders 
of the expansion about D = 4 dimensions. Numerical methods based on the BT theorem 
have also been developed for two- loop self-energy and vertex integrals [271 [23 [291 EO] , and 
found several applications [311 E21 E21 El]- While these applications involved two- loop 
integrals with full mass dependence, the phase-space integrations were trivial. We adopt 
the BT method in our calculation and compare it to the conventional analytic approach 
in order to investigate its performance in the computation of cross sections including 
nontrivial phase-space integrations at one loop. 
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Figure 2: Representative one-loop diagrams of the SUSY-QCD corrections to the process 
e + e~ — > tbH~ . The inner lines represent gluinos (g) as well as the various top (t x ) and 
bottom (b x ) squarks with x = 1,2 labeling the mass eigenstates. 

2 The calculation 

We consider the process e + e~ —> tbH~ and its charge-conjugate counter part e + e~ — > 
biH + . The corresponding LO Feynman diagrams are shown in Fig. [TJ Since the cross 
sections for both processes are identical due to CP invariance, we only show results for 
the tbH~ final state in the following. Furthermore, we only consider the kinematical 
regime where neither the intermediate Higgs bosons nor the top quark can be on shell 
and, therefore, we do not have to include the respective finite widths. 

The SUSY-QCD corrections arise from loop contributions containing gluinos, stops, 
and sbottoms in the propagators; see Fig. [2] for some representative diagrams. This 
virtual contribution may be classified into two-, three-, and four-point integrals. Since 
the particles in the loops are massive, no infrared or collinear singularities arise and no 
real-emission contribution has to be included. 

The calculation is performed with the help of the program packages FeynArts [35] 
and FormCalc [36]. The analytical expressions from FormCalc are post-processed and 
translated to a C++ code. The one-loop integrals are evaluated using the BT method as 
described in detail below. We perform a second calculation using analytical loop integrals 
as implemented by LoopTools [35]. This calculation is also based on FeynArts and 
FormCalc. The phase-space integration is performed with the Monte Carlo integration 
routine Vegas [37] in both cases. 

The renormalization of the strong-coupling constant a s is performed in the modified 
minimal-subtraction (MS) scheme of dimensional regularization with the squarks and 
gluinos decoupled from the running. The top quark only contributes to the running of a s 
above the scale m t , by increasing the number of active quark flavors from nj = 5 to nf = 6. 
The quark masses and wave functions are renormalized on shell, including the top-quark 
mass in the Yukawa coupling. For the bottom-quark mass in the Yukawa coupling, we 
use the running QCD MS version mfe(/x), with /i being the renormalization scale, and 
optionally perform a resummation of large SUSY corrections as described below. We use 
the running bottom-quark mass, since the pure QCD corrections contain large logarithms 
of the type log(/i/mfc) originating from the Yukawa interaction. These can be resummed 
by using the running bottom-quark mass rn~b(fi) [38| [39] . The calculation of the purely 
gluonic QCD corrections showed that the bulk of the QCD corrections may be absorbed 
by using the running bottom-quark mass [11]. Since the squark and gluino masses enter 
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only at NLO we do not need to renormalize them. 

In the MSSM, two Higgs-boson doublets, denoted by H u = H®) and Ha = 

(Hj,Hj), are required. The H® field is responsible for the generation of the up-type- 
fermion masses and the field for the down-type-fermion masses, i.e. the bottom quark 
couples to H% but not to H®. Nevertheless, the coupling of the bottom quark to the H® 
field is dynamically generated via loops [10] ■ Although this coupling is loop suppressed, 
once the H^ d fields acquire their vacuum expectation values v u /d, a large value of v u 
may compensate a small loop contribution, i.e. these effects may be considerable for large 
values of tan/3 = v u /vd- These large tan /3-enhanced contributions may be resummed to 
all orders [10] by replacing the bottom-quark mass in the Yukawa coupling as 



m b (ii) ( Am b v 



where 



2a 8 (n) ~,n T( \ 
& m b = — mgyUtan/3 I(rrig,m~ b ,m~ b ), 

6 ' C) = (a 2 -b 2 )(b 2 -c 2 )(a 2 -c 2 ) bg F + ^ bg h + C2fl2 bg ■ (4) 

Here, jx is the Higgs-Higgsino mass parameter of the superpotential. In order to prevent 
double counting, an extra counterterm of the form 



5m Yuk = mb (^Am b ^1 + — (5) 

for the Yukawa coupling is needed at NLO. The resummation formalism was extended to 
also include the dominant terms in the trilinear coupling A b j4T] . Since these contributions 
are small for our parameter values, we do not include them in the resummation. 

2.1 Numerical evaluation of loop integrals 

Within dimensional regularization, any scalar one-loop integral can be expressed as an 
integral over Feynman parameters as 



1 TT & I 



[q 2 - ml] [(q + pi) 2 - m 2 ] ■ ■ ■ [(q + p N -i) 2 - m 2 N ] 
= (4vr/i 2 ) e T{N - 2 + e)(-l) N J d5 iV _i^(a; i )" (JV " 2+e) , 

where 

•1 PX\ r x n-l 



(6) 



dS n = / dx 1 / dx 2 ■ I dx n 
Jo Jo Jo 

and V is a quadratic form in the N — 1 Feynman parameters Xj, 

V(x) = x T Hx + 2K T x + L — id. 
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The coefficients H, K, and L of V are given in terms of the momenta Pi and the masses 

Uli. 

In general, the quadratic form V can vanish within the integration region, although, 
strictly speaking, the zeros are shifted into the complex plane by the infinitesimal imagi- 
nary part iS. Since the limit 5 — > has to be taken in the end, the form given above is in 
general not suited for direct numerical integration. 

Instead, the integral can be rewritten using the BT theorem [26J before attempting a 
numerical evaluation. Applied to the case of one- loop integrals, this theorem states that, 
for any quadratic form V(x) raised to any real power (3, we have 



(x - X)jdj 
2(1 + 0) 



V 1+ ^x i ) = B-V^(x i ), (7) 



where X = —K T H 1 , B = L — K T H l K, and <9j = d/dxi. Inserting this relation into a 
Feynman-parameter integral and integrating by parts, one obtains 



/ dSnV ^ = 2B(1 + (3) (2 + n + W) J dSnV ^ ~ J d5 - 1 P q XtV ^ 



where \i — Xi — X i+1 , with X Q = 1 and X n+1 = 0, and 



V(l,xi, . . . ,x n -i) for i = 0, 

V i (x 1 , . . . , x n -{) = { V(x h . . . , x^ x h . . . , x n -i) for < % < n, 
V{x\, . . . , x n _i, 0) for i = n. 

Applied to the one- loop integral of Eq. (jH]), the first term inside the brackets in Eq. (jSJ) 
corresponds to the Appoint integral in D + 2 dimensions, while the last term is a sum 
over (N — l)-point integrals in D dimensions obtained by removing one propagator. 

Recursive application of Eq. (|SJ) allows us to express any scalar one-loop integral as a 
linear combination of terms of the form j dS^ V(xi) m ~ e with any integer m > 0. A Taylor 
expansion up to 0(e a ) then results in terms of the form J dSk V m ■ log 1+a V. For m = 0, 
the integrand still contains an integrable (logarithmic) singularity, while it is smooth for 
m > 0. Although larger values of m lead to smoother integrands, the expressions also 
grow larger due to the repeated application of the BT identity (jHJ). The optimal choice 
for m depends on the chosen numerical integration routine and its ability to deal with 
integrable singularities. 

The parametric representation of tensor integrals contains Feynman parameters in the 
numerator. The procedure outlined above can also be applied in this case, so that no 
separate reduction to scalar integrals is needed. Furthermore, no inverse Gram determi- 
nants are introduced using this approach, making the latter numerically reliable also for 
exceptional kinematic configurations. 

The BT identity (jHD still contains a potentially small factor of B in the denominator. 
Although the zeros of B correspond to the leading Landau singularities, the singular 
behavior in the vicinity of this singularity is overestimated by the factor 1/B. This may 
result in numerical cancellations in the numerator leading to instabilities. We, therefore, 
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use an alternative BT-like relation for the three-point function |42], namely 



(x - X)jdj 
2 



-i n=oo 




Q(x) 



) 



V~ l - e {x) = B- e 1 + 



E 



n\ Q(x) 



In 



n 



1 + 



B 



J n=l 



where Q(x) is defined by the decomposition of V(x) = Q(x) + B. For small values of B, 
the 1/B behavior is reduced to lni?, which is in agreement with the singular behavior 
near the leading Landau singularity of a triangle diagram. 

We use an implementation of the BT method in Mathematica and C++ [43J that 
allows for the calculation of all the appearing one-loop tensor coefficients. The Feynman- 
parameter integration is performed with a deterministic integration routine, with the 
number of integrand evaluations limited to 10 5 . We verify the results for the virtual 
corrections using the BT approach against an analytic evaluation of the integrals with 
LoopTools [36] for single points in phase-space and find good agreement within numerical 
integration errors. 

We could now perform the phase-space integration of the virtual corrections using 
Monte Carlo techniques. However, using the BT approach for the evaluation of the loop 
integrals requires a separate numerical integration for each integral at every phase-space 
point. This turns out to be very inefficient compared to the analytic evaluation of the 
loop integrals. 

Therefore, we combine the phase-space integration and the integration over the Feyn- 
man parameters into a single integration that is performed using the adaptive Monte Carlo 
integration program Vegas [37J. The adaptivity of Vegas optimizes the phase-space and 
Feynman-parameter integrations at the same time, leading to a significant improvement 
of the efficiency. All results for the BT approach shown below are obtained using this 
combined integration. 



Before we present numerical results, we fix the input parameters. We use the following 
SM parameters: 



The running MS bottom-quark mass rnj,(rnb) is used as input and corresponds to an on- 
shell mass of m£ s = 4.58 GeV. As mentioned above, we use the on-shell bottom-quark 
mass everywhere, except for the tbH ± Yukawa coupling, for which we always use the 
running MS version. Both m&(/x) and a s (/i) are evaluated at the renormalization scale 
fi = (mji± + m t + m&)/3. The relevant scale for the strong coupling appearing in Amf, in 
Eq. PJ is //& = {m g + + )/3. However, the two- loop result of Refs. [HI H3, HB] for 
Arrib develops a maximum at about /i&/3 to /U b /4, which is, therefore, a more appropriate 
scale choice in Amj. For our choice of Higgs-boson and squark masses, fi accidentally 
falls into this region. For simplicity, we therefore use fib = // in our calculation. We note 



3 Results 



a = 1/137.0359998, 
m w = 80.398 GeV, 
mt = 173.1 GeV, 




a s {m z ) = 0.1184, 

m z = 91.1876 GeV, 



(9) 
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m 


m 1/2 


A 


tan /3 


sign \i 


m H ± 


Ulg 




m b2 


m tl m t2 


SPSla 


100 


250 


-100 


10 


+ 


404.5 


607.7 


514.4 


543.8 


400.7 586.8 


SPS4 


400 


300 





50 


+ 


343.3 


734.4 


617.1 


682.5 


548.8 698.7 



Table 1: Definition of the Snowmass points SPSla and SPS4- The values of the masses 
and the trilinear coupling Aq are given in units of GeV. The left pari of the table shows the 
mSUGRA parameters the S PS points correspond to (the scalar and gaugino mass param- 
eters m and mi/ 2 , the trilinear coupling A , the ratio of the Higgs vacuum expectation 
values tan/3, and the sign of the SUSY Higgs-boson mass parameter fi). The right part 
shows the masses obtained from Softsusy via renormalization group running from the 
high- energy scale to the weak scale (the charged- Higgs-boson mass rriH±, the gluino mass 
rrig, and the sbottom and stop masses m~ bl , m~ h2 , m^, m^)- 

that a shift in the renormalization scale of a s (n) in Eq. (j3J) formally creates a contribution 
beyond the order of our calculation. 

In order to fix the MSSM parameters, we adopt the scenarios Snowmass Points and 
Slopes (SPS) la and SPS4 [471 HE], which are both derived from minimal supergravity 
(mSUGRA). While SPSla is a typical mSUGRA scenario with a moderate value of tan/3, 
tan/3 = 10, the SPS4 point has a high value, tan/3 = 50. In the SPS4 scenario, we there- 
fore expect large corrections from tan /3-enhanced contributions. The input parameters for 
both scenarios are summarized in Table [TJ They are defined at a high-energy unification 
scale. We use the program Sof tsusy2 . . 17 [49J to perform the renormalization group 
evolution from the high-energy scale to the weak scale and to calculate the spectrum from 
the SPS input parameters. The resulting charged-Higgs-boson and relevant SUSY-particle 
masses are also shown in Table [TJ The latter are defined in Softsusy2 . . 17 according to 
the modified minimal-subtraction (DR) scheme of dimensional reduction. They enter our 
core calculations only at one loop, so that a change of scheme would induce a shift beyond 
NLO. For the transition to the MS scheme, such a shift would be numerically insignificant 
[50] . In order to study the dependence on tan/3 and m#±, we vary the high-scale input 
parameters and recalculate the full spectrum. While tan /3 can be directly used as an 
input parameter to the spectrum calculation, we adjust m to obtain the desired value of 
m#±. This also results in a change of the squark masses. 

We are interested in a kinematical region where the center-of-mass energy y/s is not 
sufficient to produce pairs of charged Higgs bosons, i.e. \fs < 2m H ±, but sufficient to 
produce all final-state particles on-shell, i.e. \fs > mH±+irib + mt. For a 800 GeV collider, 
this implies that the interesting region lies in the range 400 GeV < m#± < 620 GeV. 

In Fig. [31 the total cross section a(e + e~ — > tbH~) is shown for the SPSla scenario as a 
function of m#± in the range y/s/2 < m H ± < y/s—mb—rrit. The cross section falls rapidly 
with increasing distance from the pair-production threshold at m#± = y/s/2 = 400 GeV. 
Before resummation, the SUSY-QCD corrections range from —22% to —24% and exhibit 
only a weak dependence on m#± . After resummation, the residual corrections still amount 
to about —15%, which implies that the tan /3-enhanced contributions are not actually 
dominant for these values of parameters. 

Figure S] shows the m#± dependence for the SPS4 scenario with tan/3 = 50. In this 
case, the SUSY-QCD corrections are much larger, ranging from —40% to —55%. However, 
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Figure 3: Total cross section a(e + e~ — > tbH~) (left panel) and relative corrections 5 in 
percent (right panel) depending on the charged-Higgs-boson mass mn± for a center-of- 
mass energy of ^/s = 800 GeV. The other MSSM parameters are taken to have their 
SPSla values. Both the nonresummed and resummed cross sections are considered. 
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Figure 4: Total cross section a(e + e~ — > tbH~) (left panel) and relative corrections 5 in 
percent (right panel) depending on the charged Higgs-boson mass vnn± for a center-of- 
mass energy of \fs = 800 GeV. The other MSSM parameters are taken to have their 
SPS4 values. Both the nonresummed and resummed cross sections are considered. 
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Figure 5: Total cross section a(e + e~ — > tbH~) (left panel) and relative corrections 5 in 
percent (right panel) as a function of tan /3 for a center- of -mass energy of yfs = 800 GeV. 
The other MSSM parameters have their SPSla values. Both the nonresummed and re- 
summed cross sections are considered. 

the bulk of these corrections is absorbed by the resummation. The size of the residual 
corrections ranges from —8% to —17% and is, therefore, similar to the SPSla scenario. 

In Fig. [5J the tan/3 dependence of the total cross section a(e + e~ — > tbH~) is shown 
for the SPSla scenario. The cross section grows strongly with increasing value of tan/3. 
Without resummation, the relative corrections also show a strong tan /3 dependence and 
range from —24% at low value of tan /3 up to —55% at high value of tan 0. Most of these 
corrections are absorbed by resummation; the remaining corrections range from —10% to 
— 15% and show only a relatively weak dependence on tan/3. Since we change tan /3 by 
modifying the high-energy input parameters, we also need to change mo accordingly to 
keep the low-energy value of m#± at its desired value. This leads to somewhat higher 
squark masses at high values of tan/3, so that the loop corrections are suppressed there. 
Numerically, we find that the competing effects of the loop suppression due to heavier 
squark masses and of the increase in the bottom Yukawa coupling lead to a maximum 
in the relative corrections at tan /3 « 45. The difference between the nonresummed and 
resummed NLO cross sections measures the size of the higher-order contributions included 
by performing the resummation. For high values of tan /3, these contributions are of similar 
size to the residual SUSY-QCD corrections underlining the need for resummation in this 
region. 

We note that the analogous analysis of the SUSY-QCD corrections to the associated 
production of a neutral Higgs boson with a bb pair in e + e~ annihilation with and without 
resummation of tan /3-enhanced contributions yielded similar results [51]. 



9 



m H - [GeV] 



a. 



BT 



[fb] 



a 



.ana 
virt 



[fb] 



410 
450 
490 
510 
550 
590 



-5.154(7)-10~ 4 

-8.46(1)-10- 5 

-2.245(3)-10~ 5 

-1.129(2)-10~ 5 

-2.152(3)-10~ 6 

-1.409(2)-10~ 7 



-5.148(4)-10~ 4 

-8.462(6)-10- 5 

-2.245(2)-10~ 5 

-1.1295(8)-10~ 5 

-2.151(1)-10- 6 

-1.4094(9)-1Q- 7 



Table 2: Contribution of the virtual corrections to the total cross section of the process 
e + e~ — > tbH~ for a center- of -mass energy of a/s = 800 GeV using the BT and the ana- 
lytical approaches. The SUSY parameters are chosen according to the SPSla scenario. 



a[fb] 



-0.00026 
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-0.00032 
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Figure 6: Contribution a virt (e + e~ — > tbH~) of the virtual corrections (left panel) to the 
total cross section and the integration error normalized to the virtual cross section (right 
panel) as functions of the common squark mass rriq for a center- of-mass energy of y/s = 
800 GeV using the BT and the analytical methods. The remaining SUSY parameters are 
chosen according to the SPSla scenario. 



3.1 Comparison of numerical vs analytical loop integral evalu- 
ation 

Let us now compare the BT method with the conventional analytical evaluation of the 
loop integrals. In Table [2j we compare the results for the virtual cross section from 
Fig. |3] evaluated using the BT method and the analytical loop integrals. The remaining 
parameters are again chosen according to the SPSla scenario. Both methods use the 
same number of points in the numerical integration performed by Vegas, although the 
dimension of the integral is higher for the numerical approach due to the additional 
Feynman-parameter integrations. The runtimes of both approaches are also similar. A 
comparison of the results shows that the agreement between both methods is good, while 
the errors are typically a factor 1.5 to 3 larger for the BT method. Detailed comparison 
of the results from Figs. H] and [5] shows a similar behavior. 
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Evaluating single loop integrals using the BT method, we observe that the errors 
and the runtime may depend sensitively on the presence of thresholds of the particles in 
the loops. Below all thresholds, the integrand is positive definite in Feynman-parameter 
space, so that the loop integral is real and the numerical integration is very fast. Once 
a threshold is crossed, the integrand becomes indefinite, the integrals develop imaginary 
parts, and the deterministic integration routine needs many more integrand evaluations 
to reach a reasonable accuracy. In order to find out if this behavior carries over to 
the combined phase-space and Feynman-parameter integration, we calculate the virtual 
corrections keeping the center-of-mass energy fixed while varying the squark masses. In 
this way, we can study the behavior in the vicinity of the various squark-squark thresholds 
in the loop diagrams. Figure [6] shows the cross section and the integration error for the 
virtual corrections as a function of a common third-generation squark mass trig = mj = 
mi 2 = mi = rn~ b2 . For the remaining parameters, we again use the SPSla values. If a 
common squark mass is used, the highest threshold in the loop diagrams is the squark- 
squark threshold at irtq = y/s/2 = 400 GeV originating from diagrams of the type shown in 
Fig. Efa) and (e). Cutting both squark lines in diagrams (b) or (c), gives another squark- 
squark threshold at rriq = mn± /2 ~ 200 GeV. The squark-squark threshold from diagram 
(d) is at rrig = St&/2, which is constrained by kinematics as < i/i — m#± ~ 400 GeV. 
Here, s t b denotes the invariant mass square of the top-bottom system. Since this is close 
to rriH± f° r our choice of parameters, the virtual charged Higgs boson can become almost 
on-shell resulting in an enhancement of diagram (d). The corresponding threshold below 
rriq ~ 200 GeV may appear to be dominating the behavior of the total cross section in 
Fig. El 

The error in the BT method shows a similar qualitative behavior to the error in the 
analytical method, although it is a factor of 2 to 3 larger. In particular, no increase 
of the error below the threshold at rriq = 400 GeV can be seen, and the increase in the 
error below 200 GeV is only moderate. The strong variation of the errors across thresholds 
observed when calculating single one-loop integrals in the BT method using a deterministic 
integration routine is almost completely absent when using the combined phase-space and 
Feynman-parameter integration. 

4 Conclusion 

The discovery of charged Higgs bosons would be direct evidence of physics beyond the 
SM. Since charged Higgs bosons are essential ingredients of all SUSY models and, more- 
over, future colliders have the potential to detect them, it is worthwhile to study their 
phenomenology. If the center-of-mass energy of the collider is not sufficient to produce 
charged-Higgs-boson pairs, the process e + e~ — > tbH~ + c.c. is of particular interest. We 
calculated the cross section of this process in the MSSM to NLO with regard to the strong 
interaction, by complementing the known SM-QCD corrections [ 1 1 J by the SUSY-QCD 
ones. 

As for the loop integrations, we deviated from the familiar analytical approach, in 
which these loop-integrals are expressed in terms of logarithms and dilogarithms. In- 
stead, we adopted the numerical approach of Ref. [23] based on the BT theorem and 
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demonstrated its feasibility for a multileg one-loop calculation such as ours. We com- 
pared this alternative method to the conventional analytical approach and found agree- 
ment within the numerical errors. In contrast to previous applications of the BT method 
[3T j 1321 133| 134] . we were tackling here for the first time a nontrivial phase-space integration 
and demonstrated how it can be performed simultaneously with the BT integrations. 

Concerning the size of the corrections, we found significant SUSY-QCD corrections 
ranging from —10% to —60%, which are of similar magnitude to the purely gluonic QCD 
corrections [TTJ. The bulk of SUSY-QCD corrections stem from tan /3-enhanced contri- 
butions and can be absorbed by performing a resummation of the bottom-quark Yukawa 
coupling. The residual corrections after resummation are still sizeable, ranging from —10% 
to -15%. 
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